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ABSTRACT 

Substantial advances have been made over the past decade in the 
prediction of turbulent flows. There has been extensive work in the 
development of turbulence models, particularly for use in boundary 
layer calculations. This review covers, in a common notation, the basic 
aspects of several important methods based on partial differential equa- 
tions for tne mean velocity field ar.d turbulence quantities, includinq 
the relationship between the methods and suqqestions for future develop- 
ment. New work on three-dimensional time-dependent larqe eddy simulations 
’s discussed. The emphasis is on the hydrodynamics of incompressible 
flows, but sources for consideration of heat transfer and compressibility 
are mentioned. 


★ 

This article will appear in Vol me 8 of the Annual Review of Fluid Mechanics. 


1. INTRODUCTION 


The computation of turbulent flows has been a problem of major concern 
since the time of Osborne Reynolds. Until the advent of the hiqh-speed 
computers, the range of turbulent flow problems that could be handled was 
very limited. The advances during this period were made primarily in the 
laboratory, where basic insights into the general nature of turbulent flows 
were developed, and where the behaviors of selected families of turbulent 
flows were studied systematically. For the engineer th^re were only a limited 
number of useful tools such as boundary layer prediction iK»thods based on the 
momentum integral equation witn a high empirical content. Features such as 
sudden changes in boundary conditions, separation, or recirculation could 
not be predicted by these early methods with any deqree of reliability. Very 
specific empirical work remained an essential ingredient of any engineer's 
analysis. 

Midway through this century computers began to have a major impact. 

First it became possible to handle more difficult boundary layers by complex 
integral analyses involving several first order ordinary differential equations. 
By the mid 1 960 ' s there were several workers actively developing turbulent 
flow computation schemes based on the governing partial differential equations, 
(pde's). The first such methods used only the equations for the mean motions 
but second generation methods began to incorporate turbulence pde's. 

In 1968 Stanford hosted a specialists conference designed to assess the 
accuracy of the then current turbulent boundary layer prediction methods 
(Kline, et. al 1968). The main impact of this conference was to legitimize 
pde methods, which proved to be more accurate and more general than the best 
integral methods. 

Vigorous development of more complex and supposedly more general pde 
turbulence models followed. Methods were first developed in which a pde 
for the turbulence energy was solved in conjunction with the pde's for the 
mean motion. Then, in an effort to reduce the empiricism required, models 
incorporating a pde ^elating to the turbulence length scales were studied. 



More recently there has been intense development of models involving pde's 
for all of the nonzero components of the turbulent stress tensor. 

The ability of these more complex models to produce predictions for 
the detal’ed features of turbulent flows has outstripped the available store- 
house of data against which these predictions can be compared; moreover the 
output of these programs now Includes quantities that are difficult, if not 
impossible to measure. At the same time these rapid developments were being 
made In computation, some totally new approaches to turbulence experiments 
were introduced (Laufer 1975). These centered on the observation that 
turbulent shear flows possess a remarkable degree of organization of their 
large-scale motions. New "selective sampling" techniques were introduced to 
study these structures, and a great deal has been learned. As yet the pde 
models have not made any use of the new experimental data, perhaps because 
large scale transport is not really consistent with the "local" ideas used 
in pde models. 

One new approach that appears promising, and is just beginning to be 
carefully explored, is the idea of using a very fast, very larqe computer 
to solve three-dimensional time-dependent pde models for the larqe scale 
turbulence. These would incorporate a simple model of the small scale tur- 
bulence in some < emi -empi rical way. At present these methods are in their 
infancy, but alreidy they have begun to shed some liqht on the simpler pde 
models, in some cases producing numerical values for constants used in the 
"simpler" two-dimensional steady pde models. As experience with this approach 
grows, and as machines improve, it seems quite likely that this type of 
calculation will eventually be useful at the engineering level. 

This review will outline the essential ingredients and effectiveness cf 
several levels of turbulent flow pde models: 

1. zero-equation models - models using only the pde for the mean 
velocity field, and no turbulence pde's. 

2. one-equation models - models involving an additional pde relatinq 
to the turbulence velocity scale. 

3. two-equation models - models incorporatinq an additional pde 
related to a turbulence length scale. 

4. str ess-equation models - models involving pde's for all components 
of the turbulent stress tensor. 
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5. la rge-ecdy simulations - computations of the three-dimensional 
time-dependent larqe eddy structure and a low-level rodel for the 
small-scale turbulence. 

Zero-equation models are common practice in the more sophisticated 
engineering Industries, and one-equation models find use there on occasion. 
Two-equation models, currently popular among academics, have not been used 
extensively for engineering applications, probably because one can do as 
well If not better In most problems with simpler methods. Stress-equation 
modeling Is now under Intensive development; It Is essential for handling 
the more difficult flows, and will probably become standard practice in 
industry In ten years time. Large eddy simulations are just in their infancy, 
and are serving mainly to help assess the lower level models. However, in 
the long term, large-eddy simulation may be the only way to accurately deal 
with the difficult flows that stress-equation models are presently tryina 
to handle. 

Four other reviews have appeared recently covering selected aspects 
of the subject. Reynolds (1974), in a publication long delayed in press, 
outlined the state of affairs in 1970. Vellor and Herring (1973) provided 
an overview of one-equation, two-equation, and stress-equation modeling as 
of mid 1972. Cebeci and Smith (1974) have an entire book on the subject, 
concentrating primarily on their own zero-equation approach. Bradshaw (1972) 
wrote an incisive and del i qhtful review of the interplay between model develop- 
ment and experimentation that should be mandatory rea-^nq for all students of 
the f’eld. 

The present review will concentrate on the hydr amic modeling of 
incompressible flows, but sources of insight for extension to compressibility 
and heat transfer will be mentioned. 


2. Z ERO-EQUATION MODELS 


The equations describing the mean velocity field in incompressible 
turbulent flow are well known (Tennekes and Lumley 1972); they follow from 
the Navier-Stokes equation by the usual decomposition of the velocity field 
into mean and fluctuating components, u^ = + u| , and may be written as 


U.U, 

J 1.J 


’ P P, 1 


(2uS 


1j 


' P ij K j 


(2.1a) 


U 


i ,i 


= 0 


(2.1b) 
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Here we use the cartesian tensor summation convention, in which repeated 
indices are to be summed over all three coordinates. Subscripts after 
commas denote partial di fferentiation, e.q. , and the over- 

dot denotes a partial derivative with respect to time. R^. * ujuj (-pR^ 
is the Reynolds stress tensor), and * ^(U^ . ♦ U j < ) is the strain- 
rate tensor; v is the kinematic viscosity, p Is the pressure, and 
p the mass density. Note that ■ 0 by (2.1b). 

To close eqs. (2. 1 ), additional equations mus* be provided for R. . . 

In the simplest models R^ Is described by a Newtonian constitutive 
equation of the form 

R 1j * 3 q? A 1j ' 2v T S 1j * 2 - 2) 

2 

where q * R.p and Vy Is a turbulent or eddy viscosity which must be 
prescribed in some suitable manner. The q 2 term can be absorbed in with 
p , and so need not be calculated explicitly. 

In a zero-equation model Vy is related directly to the mean velocity 
field U. . For free shear flows (jets and wakes) one makes the usual 
poundary layer assumptions to simplify (2.1). Remarkable success is obtained 
with simple assumptions of the form 

Vy = K/.Ub (2.3) 

where AU Is some appropriate velocity difference associated with the flow 
(e.g., the difference between jet centerline velocity and the velocity of 
the external flow), and b is a length scale chara .terizinq the width of 
the jet. The constant K may vary from flow to flow, but is typically of 
the order 0.05-0.1. In this model the turbulent viscosity is constant acioss 
the shear layer at any given downstream station (see Schlichtinq 1968.) A 
similar sort of assumption also works very well in the outer (wake) region of 
turbulent boundary layers. 

In the wall region of a turbulent boundary layer is is essential to 
consider the cross-stream variation of the turbulent viscosity. Outside of 
the viscous region a commonly used form is 

v T * xu*y (2.4) 
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Here < Is the "Karman constant," (approximatel ** 0.4), u* Is the "shear 
velocity," u„ ■ vt /p where i Is the wall shear stress and y ■ x„ Is 
the distance from the wall. Very close to the wall, where viscous effects 
are Important, success has been had with simple modifications of (2.4) 
that reflect the effect of the wall In suppressing turbulent transport, 
for example 

Vy * *u*y(l - e* y /A ) 2 (2.5) 

where y + * yu*/v , and A f Is an empirical constant. 

Alternatively, many have used the "mlxlnq lenqth model," which can be 
generalized by 

V T ' <*■«> 

where l Is the "mixing length." In the wall region of a turbulent boundary 
layer, but outside of the viscous region, the velocity field Is known to 
behave as 


3M = 

3y <y 


(2.7) 


where U = U-j Is the flow velocity parallel to the wall. 

This is the only important element of U. With i. * i n the v/a 1 1 region, 

(2.4) and (2.7) are eguivalent. 

Patankar and Spalding (1970) were among the first to document boundary 
layer computation methods of this type, and now make programs available on 
a commercial basis. More recently Cebeci and Smith (1974) devoted an entire 
book to the subject, emphasizing their own particular computational models 
and processes of this general type. A Stanford group under W. M. Kays and 
R. J. Moffat has been working with these methods for several years, with the 
distinct advantage of doing this in parallel with their comnrehensi ve ex- 
perimental program on turbulent boundary layers with wall suction, blowing, 
pressure gradient, and heat transfer. Their own particular model is 
certainly one of the most advanced of this type, and I have cho* -0 " to ielve 
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into it in more detail to illustrate the empiricism and capabilities of 
such methods. Their present program is called STAN-5, and is available 
upon request for reproduction costs (Crawford and Kays, 1975). 

The boundary layer simplifications of (2.1) reduce 


ax >y 





(u ♦ 


1 Ml 

■A 


(?. 8 ) 


2 

where we have used s (U,V,0), * (x,y,z), and p # * p/o ♦ q /3 . In a 

boundary layer calculation, p*(x) is derived fro..i the pressure distribution 
applied by the external flow. STAN-5 uses (2.6) specialized to boundary layer 
flows , 

\> T = t 2 


au 


(2.9) 


In the outer region it uses ? = >A 99 ’ where A 99 the thickness 

of the boundary layer to the point where U is 99’’' of the free-stream 
velocity U^. The factor > is provided with a dependence on the momentum 
thickness Reynolds number R, = 0U (u /u in order to better predict low 
Reynolds number flows, 


\ 


mi n 


0.085 

0.25R,, -0,25 (1 - 67.5 F) 


( 2 . 10 ) 


Here F is a wall layer blowing parameter, V o /U ij , where is the 

velocity of injection into the flow throuqh the wall. 

The inner regions are handled by assuming that 
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ry(l 



( 2 . 11 ) 


with y - 0.41. The parameter A + is given as a complicated function of 
both the pressure gradient and blowing rate, shown in Fiq. 1. There 
v o » V Q /u. , and p + B (dp/dx)(v/pu*^) . An empirical fit to Fig 1 is used in 
STAN-5. The parameter A + determines the thickness of the viscous region, 
this will not change suddenly if p + or v* changes suddenly; to accomodate 
this delay, STAN-5 uses a "lag" equation. 
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where Is determined from Fig. 1, and x* ■ xu # /\ . In handllnq the 

heat transfer problem, sir" ar models and empiricism are required; for jetalls 
sre Crawford and Kays (1975). 

For a particular flow of Interest, U,(x) and p(x) are known, and a 
"starting" profile U(x o ,y) must be prescribed. The numerics are actually 
executed In STAN-5 using the stream function as a dependent variable and the 
mean vorticity as Independent variable, as In Patankar and Spaldlnq (1970). 

The mesh points are closely spaced In the wall reqion, and then expand out 
away from the wall. 

The resulting velocity distributions, temperature distributions, skin 
friction, and heat transfer are typically in excellent agreement with experi- 
ments, except for layers very close to separation. Figure 2 shows one of 
the greater triumphs of the STAN-5 moc'»l, the heat transfer predictions for 
a turbulent boundary layer subjected at first to strong blowing, which Is 
removed midway through a section of very strong acceleration, which In turn 
is terminated downstream. The rapid changes in heat transfer coefficient that 
accompany the cessation of blowing and acceleration are extremely diffirult 
to predict; every element of the empiricism reflected above Is essential to 
the success of this calculation. 

Two other groups are experienced in the use of zero-equation methods 
for a wide variety of problems. The first is that of T. Cebeci and A. M. 0. 
Smith at the Douglas Aircraft Corporation. They have extended their calcula- 
tions to compressible flows, flows over axisymmetrlc bodies and bodies with 
longitudinal curvature, and have done extensive calculations on aircraft 
wing and body systems. Their particular model, as well as their numerical 
technique, is outlined in detail in their book (1974), which is highly 
recommended to potential user of zero-equation methods. Cebeci (19/5) 
extended their procedures to three-dimensional turbulent boundary layers. 

A second group is that at Imperial College, under D. B. Spalding. Patankar 
and Spalding's book (1970) describes their zero-equation approach for turbuler 
boundary layers, and another book by Gosman et. al (1969) describes their 


node 1 1 nq of rec Irculating flows. The most finely tuned zero-equation model 
for boundary layers Is probably the STAN-5 program developed at Stanford as 
ar. extension of the Patankar-Spaldino approach (Crawford and Kays 1975) 
Zero-equation models like STAN-5 are extremely useful In enqlneerlnq 
analysis. However, they fall to handle some Important effects, such as 
strong surface curvature and free-stream turbulence, all Important on 
turbine blades. Nor are they accurate near separation points, or In boundary 
layers subjected to extremely strong accelerations. The more advanced models, 
which incorporate a pde for the turbulence kinetic energy, were orlqlnally Intro- 
duced In the hope of providing additional generality and at the same time to 
reduce the extensive empiricism that Is essential to success In a zero- 
equation model . 

3 . ONE -EQUATION MODELS 

An equation describing the dynamics of the turbulence kinetic enerqy 
can be derived from the Navier-Stokes equations by simple manipulations 
(Tt .ekes and Lumley 1972), 

q 2 ♦ u jq 2 fj 3 20-e) - J j , j (3.1) 

Here (P 3 ' R ij U 1 j the rate production of turbulence enerqy. 

e 3 2vs"j|sV! is the rate of energy dissipation, and Jj 3 (u!u^u j ♦ 

^p'uj - 2vu^Tlj) is the twice diffusive f1 ux of turbulent kinetic energy, 

all per unit of mass. We use sjj 3 ^(uj ^ + uj ^ ) . 

Alternatively, (3.1) can be written with r replaced by the "Isotropic 

dissipation" 3 vuT .u! . and J. replaced bv J.* « uTuTul + ^n'u . ' - 
2 1 J 1 »J J j 1 1 j P v J 

'^q . . This second form is appeal inq because of the direct appearance of 

the gradient diffusion of q by v in J^* . Some authors have Incorrectly 

termed the dissipation. At high Reynolds numbers the isotropy of the 

small scale turbulence renders ^.= c , but this is not true at low Reynolds 

numbers, or near a wall. 
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In one-equation turbulence models, (3.1) form* the basis for a model 
equation 'or the turbulence velocity scale q . Typically (2.2) Is used 
as a constitutive equation, and the turbulent viscosity Is modeled by 

v T * c 2 q« (3.2) 

The length scale i Is prescribed, much as In the zero-equation approach 
typified by STAN-5. The dissipation and transport are modeled In terms of 
the scales q and t . It Is well known that, at hlqh Reynolds numbers, 
the rate of energy dissipation Is controlled by Invlscld mechanisms (non- 
linear Interactions that cascade enerqy to smaller scales) and that the 
small scale motions adjust In size to accomodate the Imposed enerqy dissipa- 
tion. Hence, by dimensional analysis 

- c i 3 /t (3.3) 

The diffusive flux is usually treated by a gradient diffusion model, 

JJ * - ( c 4 uy + y ) q 2 *j (3.4) 

STAN-5 has the capability of Incorporating this one-equation model for 
boundary layer analysis. The zero-equation approach described above Is used 
for y + ' 2A + ; for y + > 2A + (3.1 )-(3.3) are emnloyed, using (2.10) and 
(2.11) to prescribe £ . Guidance in selection of the constants is obtained 
using the well-known fact that, immediately outsiue of the viscous layer, 

(2.7) holds, and the production and dissipation terms are essentially in 
balance. Using (2.7), (2.9), and (3.2), in this reqion c? * u*/q . Setting 
(P-^ s 0 in this region, one obtains c 7 * (u*/q) * c ? . STAN-5 uses 

^ ^ 2 p 

c 2 = 0.38, = 0.055, suggested by experiments which show qVu£ ~ 7 in 

this region, and c^ a 0.59 , which was determined by comparing 

calculations with the one-equation model with those of the zero-equation model. 

2 

As a "boundary" condition on the q calculation, which is carried out only 
for y + > 2A + , STAN-5 requires that q^ be such that v, at y 4 = 2A + 
matches Vj generated by the mixinq length model (2.9) at this point. Kays 
and his co-workers have used this model to explore the effects of free-stream 
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turbulence on boundary layer heat transfer (Kearney •*. al 1970) and presently 

are using the model to st’jdy the effects of rapid chanqes In free-stream 

conditions ("non-equi 1 Ibrlum" boundary layer behavior). 

Norris and Reynolds (1975) oroposed a one-equation model t u «t shows 

promise as an alternative to the highly empirical A* correlat v and 

empirical laq equation needed If one Is to qet qood results In the viscous 

region. Their Intent was to develop a one-equation model that Is valid 

right down to the wall. Noting that at low Reynolds numbers the dissipation 

2 2 

should scale as vq /£ , they use 




1 


♦ 


_i 5 

qt/v 


(3.5) 


They argue that the length scale should do nothlnq special In the vlscojs 
region, but should behave like £ * *y right down to the wall. Near the 
wall, q y , and so (3.5) near the wall becomes^" c^vq 2 /^ and ^ 
approaches a const nt as y - 0 . This Is Indeed the proper physical 
-‘havior of the dissipation. Finally, they use (3.4), but assume that the 
turbulent transport Is suppressed by the presence of the wall, and hence 


Vy ■ c^qf(l - e’ 0 ^^') (3.6) 

Note that this produces v- ^ y 4 as y •* 0 . At the wall (3.1) becomes 
voq 2 /ay 2 * 0 , which requires CjCg/tc 2 « 1 if q - y near y * 0 . 
Having established c^ , this determines c^ . Finally, a value for 
can be estimated from the known behavior for a flat-plate boundary layer, 
and they used Cg s 0.014 . 

Norris and Reynolds applied this model to channel flow with blowlno 
from one wall and equal suction on the other. For £ they used a smooth 
fit between £ * 0.4y near the wall and £ * 0.136 in the center, where 
6 Is the channel half-width. The mean velocity profiles calculated In 
the wall region, and the change in skin friction over the no-blowing case, 
are in excellent agreement with the corresponding data for flat-plate 
boundary layers. Since the main effect on A + is that of v o + , and the 
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Norris-Reynolds model seems to handle that quite well, It does seen likely 
that it will handle the pressure gradient system as well. A boundary layer 
version of this model Is now beinq prepared to study this conjecture. 

A similar approach was adopted by the Imperial Colleqe group, reported 
by Wolfshtein (1969). However, Wolfshtein allowed the length scale to 
depart from <y In the viscous region, but kept the same behavior (3.3) 
for the dissipation. When placed in comparable forms, the constants used by 
Wolfshtein and by Norris and Reynolds are quite similar. 

Norris and Reynolds discovered an interestlnq aspect of the behavior 
of their model. They solved the channel flow equations by quessinq a 
wall dissipation, integrating outwards from the wall, and then adjusting 
the wall dissipation until the proper conditions were satisfied at the 
channel centerline. The calculation proved enormously sensitive to the 
wall dissipation, and a double precision integrating scheme had to be used. 

O 

The guessed dissipation had to be within one part in 10 of the proper value 

before tne calculation could even continue to the centerline (if the value 

2 

was further off, q either blew up quickly or went neqative). This very 
narrow window meant the! a wide variety of centerline conditions could be 
satisfied with almost identical distributions of mean velocity and kinetic 
energy in the viscous regions; computationally the model confirmed the 
concept of the law of the wall! 

Most workers have abandoned one-equation models in favor of two-equation 
or even stress equation models. However, it may be that one can do better 
with this sort of one-equation model in most flows of interest, fo- it may 
be easier to specify the length scale distribution than to com- 
pute it with a pde. This would be particularly true if the lfnqth scale 
really should be governed by the global features of the flow through an 
i ntegro-di fferential equation. Hence, further study of extended one-eguations 
models is encouraged. 

Mell^r and Herring (1973) discuss some of the earlier work on one 
equation models, citing numerous references of particular calculations. 

The serious student of this subject will find their review particularly 
useful as a resource for computational examples. 
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4 . TWO-EQUATION MODELS 

In attempts to eliminate the need for specifying the turbulence lenqth 
scale l as a function of position throughout the flow, several workers 
have explored the use of a second turbulence pde which In effect gives l . 
The groups at Imperial College and at Stanford both experlmerted with ad-hoc 
transport equations for l , with no real success. However, success has been 
had by both groups and others uslnq a model equation based on tne exact equa- 
tion for the Isotropic dissipation^- ; this equation can be developed from 
the Navler-Stokes equations by appropriatt differentiation, multiplication, 
and averaging, and Is 

• 

& tu j^j * -"-Vi (4 -' a) 


here 


‘ 2vu i * 2 u 2 u i,jj u i,kk 

* 2v ( u i.j u l,k U j,k * u i,k u ).k u 1 ,j) * 2 '=Kk\,k < 4 -" >) 

■ vu i ,k u t ,k u 3 * ^"j.kP'.k - (4 - ,c) 


H. represents the diffusive flux of In the j direction. 

J 

The systematic workers have insisted that their two-equation models 
first describe properly the decay of isotropic turbulence, and then have 
worried about the behavior of their models in homoqeneous flows where the 
transport terms vanish. For the isotropic decay problem, (3.1 ) and (4.1) 
reduce to 


q 2 = - 2 = -W (4.?a,b) 

W is a scalar for wMch a closure assumption is needed. In this problem 
W must be a function of the only other variables around, and , and 
from dimensional arguments must be (at hi qh Reynolds number) 

W » c ? £:Vq 2 (4.3) 
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The exact solution for the decay is 


q 2 » q£ (ln/a)’ n ^ (ln/a)’^ n+1 ^ (4.4a,b) 

a « nq 2 /(2^ ( ) n * 2/(c 7 *2) (4.4 c ,d) 

Mere and ^ are the initial values. Early experiments sugaested 
n * 1 , which gives c 7 * 4 . Comte-Bellot and Corrsln (1966, hereafter 
denoted by C-BC) took special care to obtain better Isotropy, and their 
data reveal n values in the range 1.1-1. 3. Lumley and Khajeh-Nouri (1974b, 
hereafter denoted by LK-M2 ) suqgested that sliqht anisotropies are responsible 
for these differences, and proposed a righer-order model to take this into 
account. But this theory does not explain the different values observed in 
truly Isotropic decay, as revealed in Table 3 of C-BC. It seems more reasonable 
that the structure of the low wavenumber portion of the soectrum is responsible 
for these differences. 

The influence of the low wavenumber spectrum on n can be shown usinq 

the spectrum of Fig. 3, following a similar analysis of C-BC. The low wave- 

number part of the spectrum is assumed to be permanent, and the hioh wavenumber 

portion moves down as becomes smaller. The peak, which corresponds to the 

energy-containing scale, occurs at wavenumber k. . To the left of the peak 
m 4 

we take E s Ak ; it is known that E ^ k for k -► 0 , but this might not 
include the energy containing range and so we allow a less gradual growth in 

4 

this range. For k < k Q the k behavior would obtain, but we shall not 

need to deal with this region. To the right of k. we use the Kolmogoroff 

-5/3 L 

inertial sub-range spectrum E o, k ' . The constant a is universal for 

this spectrum, and has a value of about 1.5. In the inertial sub-range energy 

is transported up the wavenumber scale by non-lineat interactions, and the 

spectrum is controlled solely by the rate at which energy is being processed 

upscale (i.e., by the dissipation ^ ) . At hiqh wavenumbers viscosity is 

important, but this range does not contain significant energy and need not be 

considered here In detail. It is a simple matter to calculate the energy 

contained in this model spectrum from q 2 /2 * / E(k)dk, assumino k <<k. «k . . 

J o L d 

One finds 

*!)V 2/3 ^ /3 <«•*> 
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It is interesting that the form of the large-eddy spectrum enters through m , 


but its strength (A) does not. Egn. (4.5) shows that the length scale of the 
energy-containing eddies is g^/^- (compare 3.3), and hence the time sc le 
is q 2 /\>. 

Matching the two portions of the spectrum gives 

k L ■ [aW /3 /A] 3/(3m * 5) 

Then, using (4.5) 

* C [ q 2](3m+5)/(2m+2) (4 .6) 

2 

where C Is a constant. Substituting in (4.2a), and solving for q , one 
obtains (4.4a) with n ■ (2m+2)/(m+3) . So, m * 4 gives n = 10/7 , m * 2 
gives n * 6/5, and m » 1 gives n * 1 . 

It Is clear that the details of the low wavenumber portion of the spectrum 
are instrumental in determining n ; since these details are in no way 
represented by the scales q and ly , there is no way that this model can 
e>actly predict the decay of laboratory grid turbulence. However, it Is 
possible to make a fairly rational choice of c -j . We really should expect 
the model to work only when the large-scale structure is devoid of any scales, 
i.e. when the larqe-scale energy is uniformly distributed over all wave vectors . 
This occurs only when $^(k_) is the same at all k_ low wavenumbers. The 
three-dimensional energy spectrum function used above is E(k) = 2nk^ . (kj , 
and represents the energy associated with a shell of wave vector space. Hence, 
in "equipartioned" large-scale turbulence, E(k) % k .On this basis we 
recommend n = 6/5 , which gives c^ = 11/3 . This is close to the value 
used by '.K-N2 and the Imperial College v/orkers. 

When strain is applied to the flow, there is every reason to expect an 
alteration in W ; something must provide a "source" of^^, and this must 
depend in some way on the mean flow. Lumley has argued that this can not come 
from the terms in W explicitly containing the mean velocity, but must come 
from the first two terms in W (see 4.1b), which are very larqe but of 
opposite sign. Lumley feels that the alteration of W by strain should be 
modeled in terms of the anisotropy of the Reynolds stress tensor. If we 

t 4 

There is no real reason to require E(k) v k , as required by analyticity 
in k^ as k -* 0 (see Hinze 1559). The box-like grid certainly could create 
a di rectional ly-dependent dE/dk^ for k_ * 0 . 
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follow this approach, and represent the anisotropy through 


b ij * (R 1 • - q 2 6 1j /3)/q 2 (4.7) 

the.i the first scalar that can be formed from the anisotropy measure is 
2 

b ■ b.,b^j . Lumley therefore proposes 

W - (C 7 - Cgb^V ( 4 - 8 ) 

LK-N2 use c 7 r 3.73 and c Q ■ 30. 

° 2 

In a two-equation model b must be produced from the constitutive 
equation (2.2), with Vy given by 

v T ■ c g q 4 /^- (4.9) 

A 

To match (3.2) and (3.3), c Q 3 c 7 c 7 . Then, b.. » 2c Q q^S,./i^, 

2 242jcx y) t s 1 J y 1 J 

fc 3 4c g q S , where s \ s u s u • The turbulence production is 

^ 3 2c g q 4 S 2 /^. , and hence b 2 3 2c g Hence, in this model (4.3) 

may be written as 

W - (c 7 - c 1Q $/b.) &7q 2 (4.10) 

where c 1Q 3 2c g c g . 

Using Lumley's value of c g = 30 and the other constants qiven earlier, 
c^q = 1.25 . The qroup under B. E. Launder at Imperial College have explored 
two equation models extensively, using forms equivalent to (4.10) with c^ = 3.1 
It seems most desirable to determine Cy q by reference to experiments 
in nearly homogeneous flow, where the transport would not confuse the isswp. 
There are two types of such flows, those involving pure strain and those 
involving pure shear. Tucker and Reynolds (1968, hereafter denoted by TR) 
and Marechal (1972) studied the pure strain case; Champaqne, Harris, and Corrsin 
(1970, hereafter refered to by CHC) and Rose (1966) studied homoaeneous 
shearing flows. In 1970 (4.10) was proposed as a generalization of models 
used by Launder and others, and the constants were evaluated by reference 
to the TR and CHC flow (see Reynolds 1974). For that evaluation c 7 3 4 was 
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used. Recently L. H. Norris and I repeated the evaluation for the preferred 
value of c 7 * 11/3 . We carefully evaluated the prod.jctlon term from the 
data In these two flows, and used this as Input to (4.10) . The q history 

was carefully differentiated to get an Initial value for ^ , the ^ and 

2 

q equations were solved simultaneously by an accurate forward-difference 

2 

Integration, and the q histories were compared with the experimental data. 

We found that c^g ■ 2 gives excellent agreement In both flows, as found in 
the earlier work. Hence, if one elects to use (4.10) In any model, the choices 
c^ * 1 1/3 , c^g ■ 2 are recommended. 

At this point we have a two-equation model that can be tested against 
the homogeneous TR and CHC flows. In a prediction the R^ and hence (? 
must be derived using the constitutive equation (2.2) with (4.9) . Remarkably 
good results are obtained for the TR flow with Cg * 0.025 . As noted below 
(4.9) , gives Cg s 0.020 using STAN-5 constants. With this value the two- 
equation model underpredicts (? In the TR flow, and does not produce enough 
anisotropy in the Reynolds stresses. When applied to the CHC flow, the two- 
equation model fails miserably in prediction of both shearing and normal stresses. 
A weakness of (2.2) is that It forces the principal axes of R^ and 
to be alighned. This is true In pure strain (the TR flow), but not 
true in any flow with mean vortlcity (e.g., CHC) . One Is tempted to try 
a modified constitutive equation (see Saffman 1974) 

R 1j “ V 6 1j " 2v T S 1j ’ c ll 4 ^ S 1k ft kJ + S jk fz k1 ^ 

where P.^ = ^ - Uj ^) Is the rotation tensor. In a two-equation mcdel 

£ could be expressed in terms of q 2 and ^ . Eqn. (4.11) does produce the 
right sort of normal stress anisotropy in shear flows, but the new terms don't 
alter the shear stress, and hence (4.11) works no better than (2.2) for the 
CHC flow. Two-equation models also fail to predict the return to Isotropy 
after the removal of strain, or the isotropizinq of grid-generated turbulence 
(C-BC). This failure arises because of the need for a constitutive equation 
for the R^j . Thus, one should not really expect two-equation models to be 
very general, although they might be made to work well with specific constants 
in specific cases, such as boundary layers. 
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In spite of these difficulties with models based on constitutive equa- 
tions, their simplicity makes them attractive. Two-equation models have 
been studied by a number of groups, and It Is significant that these workers 
Inevitably turn to stress equation models because of the difficulties out- 
lined above. Stress equation models have their own problems, and so there 
probably Is still considerable room for development of two-equation models. 

Of particular Interest Is turbulent boundary layer separation, where anisotropy 
of the normal stresses Is known to be Important. Since (2.2) won't qive this 
properly In a shear layer, but (4.11) can, the use of (4.11) In conjunction 
with two-equation models should be explored further. 

To use the two-equation model outlined above In an Inhomoneneous flow, 
one needs to assess (or neglect) the effects of Inhomoqenel ty on W , and 
also to model the transport term Hj . Jones and Launder (1972, 1973, 
hereafter referred to by JL1 and JL2) assume that W Is not modified by 
inhomogeneity and use a gradient diffusion model for Hj , 

Hj ■ -(v ♦ c^y) &»j (4.12) 

with c^ = 0.77 . Lumley (see Lumley and Khajeh-Nouri 1974a, hereafter 
denoted by LK-N1 ) argues on formal grounds that the diffusive flux of dissipa- 
tion should depend as well on the gradients In turbulence energy, and vice 
versa. In the manner of coupled flows such as thermoelectricity and thermo- 
diffusion studied by the methods of Irreversible thermodynamics. If this is 
true, one really should use models of the form 


■ A ll q ’j ’ A 12^*j 

(4.13a) 

' A 21 q ’j ' A 22^ ’j 

(4.13b) 


Lumley and his coworkers have done this In their stress-equation model Inq, 
but as yet no users of two-equation models have adopted this approach. 

Eqn. (4.13) would allow for up-gradient diffusion of turbulence energy, a 
real phenomena in the central region of a wake, while the simpler uncoupled 
models do not. This is an area worthy of further experimentation within the 
structure of two-equation models. 
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The ^ equation model described above works fairly well at high 
Reynolds numbers, but falls near a wall where viscous effects are Important. 

JL proposed ad-hoc low Reynolds number modifications that seem to work reason- 
ably well In the wall region, and Hanjallc and Launder (1974, hereafter denoted 
by HL) proposed further modifications of theS^r equation for use with their 

stress-equation model. Clearly the W term has to be modified, for in the 

2 -5/2 

"final period" of decay of Isotropic turbulence q ^ t Instead of 
t"^ . If the turbulence Reynolds number Ry * q 4 /(fcv> Is small , then the 
inertial terms In the dynamical equation are unimportant, and In Isotropic 
turbulence W Is dominated by the second term In (4.1b). At low Ry , 
fcr -v vq 2 /l 2 , so l ^ (vq 2 /J^) 1/2 , and W ^ v 2 q 2 /fc 4 ^ & 2 /q 2 • Hence, at 
low Ry , W * Cy^/q 2 , which is of the same form as Lhe high Ry behavior 
(see 4.3). Setting n ■ 5/2 In (4.4), c* B 14/5 , which Is consistent with 
the models of HL and JL. A smooth transition between c^ and Cy Is needed; 
a form similar to that used by JL and HL but consistent with c^ * 11/3 and 
c* « 14/5, Is 


W * y ^(Rylfc-W (4.14a) 

where 

* 1 - g exp [-(Ry/12) 2 ] (4.14b) 

Remember that this Is just for the part of W that Is non-zero In homogeneous 
Isotropic turbulence. 

Eqn. (4.14) presents problems near a wall, where -*■ const, and q 2 - 0 . 
Launder and his coworkers get around this by ad-hoc modifications of their 
model equations. HL replace & 2 by WW where ^ - vfoq/ax^) 2 . 

Unfortunately they refer to as the isotropic dissipation, for some reason 
confusing it with^* . In spite of this semantic problem, their assumption 
does seem to work in boundary layers. However, this reviewer would prefer 
an approach in which - const, as y - 0 , which physically Is 
correct. 
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An alternative approach to handling this part of W near a wall Is 


11 

T 


f,(R T )(l - e-'ISiy/^V 


(4.15) 


This gives W -*■ const, as y -* 0 . The use o 1 (4.15) should be explored. 

The third and fourth terms on the right In (4.1b) vanish at hi qh Ry 
because of the sn«ll scale Isotropy. HL planned to Include these at low Ry 
by lumping them with the first two terms In (4.1b) by further modification 
of f^ . However, they found that this was not necessary. The last 
term In (4.1b) was neglected by JL. It was modeled by HL In a complex way 
Involving products of two second derivatives of the mean velocity and the 
Reynolds stresses. 

JL2 used the two-equation model to study a limited number of boundary 
layers. Including the "difficult" flow shown In Fig. 2. The predictions of 
their model are seen to be noticeably less accurate than those of the STAN-5 
one-equation model shown. 

One difficulty with using the ^ equation as the basis for a second 
model equation has escaped the model developers. This arises from the second 
term In (4.1c), the pressure gradient-velocity gradient term In the transport 
Hj . Since the pressure field depends explicitly upon the mean velocity field 
(see |5) , mean velocity gradients can explicitly give rise to transport 
This could be an extremely important effect, especially near a wall. The 
omission of this consideration would seem to be a serious deficiency in all 
equation models that have been studied to date. 

Other two-equation models have been heuri stlcal ly conceived. Of these 
the most well developed is the Saffman-Wi 1 cox (1974, hereafter denoted by SW) 
model. Instead of a ^ equation they use an equation for a "pseudovorticity 


•2 2 
a + u j 


(4.16) 


[ a Vu7^TT • ^ 

* [< v * ov t) 2* j ] . 

• J 

In conjunction with this they use the q equation (3.1) with 

9= a*^ 2 S^*q 2 /2 S*q 2 n/2 (4.17a,b) 
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2 

They use (3.4) for the q transport, setting c 4 ■ o* , and for \y they set 

v T » q^/(2£) (4.18) 

The constitutive equation (2.2) Is used to provide R^ for the mean momentum 
equations. Their recommended constants are a ■ 0.1638, a* ■ 0.3, 6 * 0.15, 

6* ■ 0.09, o ■ 0.5 , c* ■ 0.5 . 

The production term as given by (4.17a) Is Inconsistent with the R^ 
constitutive model; this seems to be an Internal Inconsistency In the model, 
but it may In fact be a strength. The 9 model Is based on the experimental 
fact that the structure of the turbulence In the wall region of a boundary 
layer Is essentially Independent of the strain rate, and hence ^ should 
be proportional to q . Hence, the SW model Is a curious blend of the 
"Newtonian" and "structural" alternatives (Reynolds 1974). 

For Isotropic turbulence decay the SW model equations may be solved 
exactly. The high-Ry behavior, q^ *v. , Is obtained If f?*/8 * 3/5 , 

as suggested by SW. I recently tested the SW model against the TR and CMC 
flows, using "starting" values for u carefully calculated from the Initial 
q^ decay rate. In neither case were the results at all Impressive, ;ioreover, 
the SW model does not display the proper decay of Isotropic turbulence at low 
P T . Therefore, it does not appear that the SW model Is or can be any more 
general than any other two-equation model. Indeed, both Saffman and Wilcox 
are Independently exploring stress equation models (Saffman 1974, Wilcox 1975), 
neither version of which presently works very well In the TR and CHC flows. 

The SW model has been tested against only a limited body of boundary 
layer flows. The model works surprisingly well In the viscous region, but 
has the troublesome point that £ must be infinity at a perfectly smooth 
wall. SW use a "large" value of £ at the wall to produce mean velocity 
curves that are in excellent agreement with experiments for smooth walls. In 
effect, SW match their solution to experimental data by judicious choice of 
the value of the wall £ . In SW they considered only zero pressure gradients 
with no transpiration. More recently Wilcox (1975) examined a few cases of 
pressure gradient and transpiration, and made a useful comparison of the SW 
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model with other two-equation models. Including JL. By jut ous selection 
of the wall value of w he could match some of the Stanford tansp' red 
boundary layer; his calculations Indicated a strong effect of blowlnq on the 
wall £ . Thus, the SW method will require a graph of the wall £ as a 
function of pressure gradient and blowing parameter, similar to F 1 q . 1. 

Wilcox also found It essential to use accurate values for the free-stream 
value of u , which he also had to carefully deduce from experimental data. 

It appears that the sensitivity of the SW model to free-stream conditions 
nay be significantly greater than that of the JL model, and certainly Is much 
greater than that of one-equation models. 

One Is led to conclude that the SW model should not be used as an 
engineering tool until such time as It has been developed much further. 
Regarding £ as a reciprocal time scale may be useful In qulding these 
developments. 


5. STRESS-EQUATION MODELS 

In turbulent shear flows, the energy Is usually first produced In one 
component and t 1, o transferred to the others by turbulent processes. 

Exact equations tor can be derived from the Navler-Stokes equations 

(Tennekes and Lumley, 1972); for an incompressible fluid, 


R lj * U k R IJ,k 


P. , ♦ T . . - D. . - J . .. . 
ij 1j 1j Ijk.k 


(5.1a) 


Here P^ Is the "production tensor," 


P 1j “ " R 1k U 1,k " R jk U 1 ,k 


8 - (R 1k 5 kj + R jk S k1 } * (R 1k\j + VW (5Jb) 


Note that P 


11 


2§* . is the "transfer tensor", 


u 


■ p p ‘ tu i.j * u j.i» 


(5.1c) 
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This "pressure-strain" tern Is responsible for energy exchange between 
components. Note that T^ ■ 0 by continuity. ; Is the "Isotropic 
dissipation tensor," 


'u 


2v 


(5. Id) 


and 


Is the diffusive flux of R^ , 


.1 . * — ( p ' u 1 6 . , ♦ n ' u ' 6 j , ) ♦ u 1 ii 1 ii ' _ >iP I X lol 

IJk P ' 1 Jk 1 J IK* 1 J k 1J ,K 


Note that J... ■ J? . 

11k k 

P^ Is explicit, but models are needed for T^., D^, and J.j k . 

In addition, one must either specify i or use a equation. We will 
first discuss the hlgh-Reyno’os number modeling of (5.1), particularly as 
applied to homogeneous flows, and then discuss the problems and status 
of extending this model to Inhomooeneous reqlons, particularly near walls 
where Ry Is smal 1 . 

Tne one fact that seems very c'ear from experiments Is that, at high 
Ry the small-scale dissipative structures are isotropic. Hence all workers 
now use 



(5.2) 


The transfer term has been the subject of most controversy and 

experimentation. In a flow without any Piean strain, this term is responsibl 
for the return to isotropy. However, in deforming flows the situation is 
much more complicated. Guidance is provided by the exact equation for the 
fluctuation pressure, derivable from the Navler-Stokes equation (see 
Tennekes and Lumley 1972), 



■ 2u i.J U j.1 




9l ♦ q 2 


(5.3) 


The source term in this Poisson equation contains two parts, each of which 
will be responsible for a part of the pressure field. The part determined 
by g^ , which involves the mean deformation explicitly, we denote by p,J , 



and the remainder by pA . Followlnq LK-N2, the explicit dependence of 
the p-j contribution to T.j can be obtained for honoqeneous fields In 
terms of the Fourier transform of the velocity field. Let 


P{ * J pfkje 1 -'•*<% 


(5.4) 


In homogeneous flows the mean qradlents are constants, so (5.3) nlves 

k 
k 


P “ 


2u t iv 


(5.5) 


Then, the part of the pressure-strain term associated with pj Is (we 
aJopt the subscript choice of LK-N2 for convenience In comparison) 

T lpq ■ Pl'Vq'Vp 1 • ff pfk) [ u f (k '>% * ^(k'lkjdk'dk (5.6) 
Uslnq (j.5) and the statistics of random transforms, (5.6) becomes 


T )pq " 2u jJ G ijpq 


(5.7a) 


where 



(5.7b) 


Eqn (5.7a) Is Identical to an expression developed by Rutta (1951) from 
slightly different arguments. 

Models for G^j pq have been proposed by Launder and Lumley and their 
coworkers. There are various const, alnts 'at G. . must satisfy. From 

i jpq 

continuity, G^j pp * 0 , G 1ipq s 0 . Also, G^^ ■ R^ q . For Isotropic 
turbulence these suffice to define G^ pq . Harjailc and Launder (1972) 
first used a model of G^ that Involved linear and quadratic terms in 
the R.j . Later Launder, Reece, and Rodl (1973, hereafter denoted by LPR) 
dropped the quadratic terms. LK-N2 also used non-linear terms, but later 
Lumley (1975) arqued that the model must be be linear in the Pevnolds stresses 
because for a field that Is the sum of two uncorrelated fields T^. . should 
be the sum of their Individual ^ . s . Lumley (1975) sought to resolve 
certain inconsistencies between the calculations and experiments by allowing 
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**1jpq t0 depend a compl Icated way on scalars developed **•<% , combinations 
of the mean deformation and b^ (see 4 . 7 ) . But this violates the condi- 
tion that the should not be changed by a sudden change In the mean 

strain-rate. If this condition Is Imposed, and we Insist on linearity In 
the Reynolds stresses, then the model (In a homogeneous field) must 

be of the form 

G IJpq ' j*TC 4 tJ 4 pq * iV * 4 tp 4 Jq * 4 1q 4 jp*j ^ 

* 1 jlV'jq * Vjp> * ¥' "jp*1q * Vl P > 

' (5.8) 

* A l [b 1J 6 pq ' ^ b 1p 6 Jq ^ b 1q 6 Jp } 

‘ ^hjp-lq * b Jq 6 1p } * b pq 6 1J ] | q2 

Using this In ( 5 . 7 a), the part of T^ explicitly related to the mean 
field must be 

T 11j " S 1j q2 " F A 1 [ R 1k S kj + R jk S k1 + I^ 4 1.l] 

4/5 7 V r ( 5 . 9 ) 

• f(t + t/i) [ R 1k n kj + R jk Q kl] 

This Is precisely the form used oy HL. 

The part of T^ associated with g^, which we denote by T^ , 
should not change Instancy when the mean deformation Is chanqed, and 
hence should not depend explicitly on the mean deformation. LK-NT Ignored 
this requirement, and allowed to depend on the rotation tensor. 

Lumley ( 1975 ) has now abandoned this position. Launder and his coworkers, 
and others, have followed Rotta In assuming 

T 21J * - A oV>U (5 -'°> 

The constant A q determines the rate of return to Isotropy. Its value 
has been the subject of much uncertainty. The TR flow Implies a value 
A q s 6 , while the C-BC data suggest a much lower value Is appropriate. 
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ML and LRR use A q • 3.0, LK-N2 use A q ■ 3.21 . LRR point. to the advantaaes 
that would be obtained If a lower value of approximately 0.6 could be used, 

In which case the behavior near a wall would be much more accurately modeled. 
Lumley and his coworkers add additional non-linear terms In the b . j , 
feeling that the rate of return to Isotropy should depend upon the degree 
of anisotropy. It does not seem that the data justify the inclusion of 
higher-order t-rms, and so (5.10) Is recommended, at least for homogeneous 
flows away from boundaries at high Ry . 

L. H. Norris and I recently studied this problem using the exact 
solution of the model equations for ' ,ie retur to Isotropy In homoqeneous 
turbulence without strain. Using (5.2) and (5.10) In (5.1), for this case 



Iqns. (4.2) aqaln describe q 2 and . The exact solution for the 
decay Is (see 4 .4) 


b ij " b ij 0 (1 + t/a)’ ( V 2)n/2 (5.12) 

where b^ o are the Initial values. Note that A q must be at least 2 

if isotropy Is to be restored. Norris and I used the data of C-BC's Table 1, 

and first simply solved (5.11) for ( A q - 2 ) . Subsequently we compared the 

solution (5.12/ to the data, using n * 6/5 . There is a great deal of 

scatter, because the anisotropies are rather small. There was absolutely 

no systematic dependence of A q on either anisotropy or Ry . Based on 

this work, we recommend A„ s 5/2 . 

o 

Kwak and Reynolds (1975) studied the TR flow ir a numerical simulation, 
and found a much slower return to isotropy than indicated in the TR experi- 
ments. However, different components return at decidedly different rates. 
Shaanan and Ferziqer (1975) carried out a similar calculation for a shear 
flow similar to that studied by CHC. In a computation the shearing can be 
removed, which can not be done experimentally. These calculations also 
showed a marked difference in the return rate for different components, 
probably because of great difference in the length scales in the Ihret 
directions. We conclude that current stress equation models will not do a 
very good job in handlinq the return to isotropy; however, the models may 
work well in flows dominated by other effects. 
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The constant A^ should be evaluated by reference to homogeneous 

flows, such as the TR and CHC flow. LK-N2 used -P.456 , which was 

obtained by a comparison with a rapid distortion analysis of homogeneous 

strain. Later Lumley (197*) argued aqalnst this approach, and settled on 

-1.23 (In a more Involved model). LRR use a value of -1.45 (their 
2 33 

c^ * -(j ♦ ^ A] )) , which they base on homogeneous experiments. I recently 
found that -1.5 is a reasonable compromise between -1 which works better 
for the TR flow and -2 which works better for CHC, and now recommend 
A 1 * -| 

Inhomogenelties greatly complicate the Tj- modeling, especially 
T^j . LRR add a complicated term Inversely proportional to the distance 
from the wall. Recently M. Acharya and I extended LK-N2's analysis for 
T^.j to a flow near a wall. We took Fourier transforms In only the x^ 
and x^ directions, and solved the ordinary differential equation for the 
transform amplitude p(y) . This leads one to a messy Integral expression 
in which T^ . j depends upon the mean velocity gradients at all points in 
the flow . In a wall region one might well expect T^j to be determined 
by a region at least as wide as the distance to the wall, and hence a complex 
integral model is really needed for such flows. This Is a very unsatisfactory 
aspect of present stress-equation modeling, and an area that should receive 
conslderabl e attention in the future. 

In addition to modifications in T^ , Inhomogeneities require modeling 


of J, 


ijk 


The qradient diffusion model is usually employed; HL and LRR set 


'ijk = - A 2]> ( R 1n R jk,n 


+ R ._Rj l + R. R . . 
jn Ik ,n kn ij,n 


(5.13) 


Hanjalic and Launder (1972) gave some justification for this form by con- 
sideration of the dynamical equation for uju 1 u^ . Lumley (1975) used 
somewhat more extensive arguments to in effect provide further justification 
for this form 
and since 


Noting that contains one pressure-veloci ty term. 


p' will have a part (pp that depends explicitly on the mean 
ve’ocity gradients, it does seem that also should be explicitly 

linear in the mean gradients, though this need has escaped notice. 
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Other modifications necessary near a wall have been sugqested by L RR . 

In particular, they propose to allow anisotropy In Dj . at low Ry , and 
have concocted a smooth transition between (5.2) and 

0 ij - 2R 1j £r/q 2 (5.14) 

which they Incorrectly imply Is exact as Ry ♦ 0 . 

Two approaches have been used In stress-equation modeling. The earlier 
work (Donaldson 1972) Involved specification of the lenqth scale and use of 
(3.3) to determined . HL used the ^ equation model outlined above in 
conjunction with the R^j equations. At this writing this work Is in a 
state of rapid development, and undoubtedly improvements will be made by 
the time this article !; released. Interested persons should follow most 
carefully the work of Launder and Lumley. It will be some time before these 
models are sufficiently well developed to be better than simpler models for 
use In engineering analysis. 

An Interesting use of stress equation models Is suqgested by a contraction 
of (5.13), 

2 

J Hk ' ‘ 2 VlK.n) (5 - ,6 > 

If this is compared with (3.4), its counterpart In the one- or two-equation 
models, an important difference is seen; Eqn. (5.14) does allow for a flux 
of q to be driven by gradients of other than q . Moreover, if the 
constitutive equation (2.2) is used with (3.2), the q flux will be driven 
by mean velocity gradients These effects are not incorporated in (3.4); 
an approach to improving the simpler one- and two-equation models miqht be 
to use the more complex stress equation model as a guide to the nature of 
new terms that should be included. 

There is a basic difficulty in this qeneral approach to turbulence 
models. One would like to model only terms that respond on time scales 
short, compared to that of the computed quantities. It is well known that 
the small scale: respond to change much faster than the large scales, and 
hence it is reasonable to express a quantity dominated by small scales, such 
as D. . , as a function of quantities dominated by large scales, such as P . . . 
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However, terms like J^ k have time scales comparable with that of R.j , 
and thus one really should not expect an equilibrium constitutive relationship 
to exist between and . In qeneral , it seems that hiqher order 

statistical quantities take longer to reach steady state than lower order 
statistics; for example, in a channel flow the "entrance length" for the 
mean velocity Is ratner shcrt, while the entrance length for the R^ • is 
known to be quite long. Any model obtained by truncation at some statistical 
order would suffer from this difficulty. What one really needs to do is 
truncate at some level of scale , and thereby take advantage of the fact 
that the smaller scales do adjust faster to local conditions. Then, by 
truncating at smaller and smaller scales, one has at least some hope of con- 
vergence, a hope that is at best dim when one truncates at hiqher and hiqher 
orders of statistical quantities that have comparable time scales. The large- 
eddy simulation described in the next section provides one avenue to a 
scale-truncation approach. 

An interesting identity that might be useful in a different approach 
to turbulence modeling is 


R. , 

ij.J 


" e 1jk + 


(5.16) 


where wj = 


e 1jk u k,j 


Stress 


is the fluctuation vorticity. When this is used in 

(2.1a), the Reynolds "stresses" disappear (except for a "Reynolds pressure" 
2 

q /?■ ), and are replaced by "Reynolds body forces" F\ = 
equation models try to model R^ , and then take their gradients. It 
might be easier to model the body forces directly. For a physical 
discussion of the F. , see Tennekes and Lumley (1972). 


6. LARGE EDDY SIMULATIONS 

This line of approach is just beginning to bear fruit. The idea is 
to do a three-dimensional time-dependent numerical computation of the large 
scale turbulence. It is impossible to compute the smallest scales in any 
real flow at high Rj (and will be forever), so they must.be modeled. 

Care must be taken to define what it is that is being computed and the 
early work was not done with sufficient care. 
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In 1973 we began a systematic program of development and exploration 
of this method, In close cooperation with NASA-Ames Laboratory. The first 
contribution was made by Leonard (1973), who clarified the need for spatial 
filtering. We now define the large-scale variables by (see Kwak and Reynolds 
1975) 

T(x) ■ /g(x-x/ ) f(x/)dx' (6.1a) 

where the fl 1 ter function Is 

G(x-x') - |^Ij-J 3 exp [-6(x-x') 2 /£ 2 ] (6.1b) 


Here a. Is the averaging scale, which need not and should not be the same 
0 

as the grid mesh width. We use this particular filter because of Its 
advantages in Fourier transformation. When this operation Is applied to 
the Navier-Stokes equation, and an expansion Is carried out, one finds 
(neglecting molecular viscosity) 


a j U 1 . J 


I F 


.1 


*2 

If (D ‘l°J ) .kk ' R 1j 


♦ 0(aJ) (6.2) 


where -pR^ are the "sub-grid scale Reynolds stresses." The unusual term 
appearing before R^ Is an additional stress-like term resulting from the 
filtering of the non-linear terms; we now call these the "Leonard terms," 
and view pA* (U^Uj) j^/24 as the "Leonard stresses." 

We have explored two models for the R i ■ . Both are based on (2.2); 
the first Is Smagorinsky 's (1963) model, 


vf 


V 


s s 

nm nm 


(6.3a) 


The second uses the rotation in place of the strain rate, 

v T = (6.3b) 

l 2 a 1 nm nm 

In these expressions and 12^ are the strain rate and rotation of 

the calculated local time-dependent large scale field. The q 2 term in 
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(2.2) Is again absorbed with the pressure. Note that the sub-grid terms 

are 0 ( ) , and hence If they are Important the Leonard stresses 

are also likely to be Important. Moreover, a difference scheme must be used 

2 

that Is accurate to 0(A) ; this Important requirement was overlooked 

d 

by many of the early workers. 

kwak and Reynolds (1974) solved the Isotropic decay problem, adjusting 
or B 2 to obtain the proper rate of energy decay. The calculations 
were started using an Isotropic field with zero skewness, but the proper 
skewness develops In only a few time steps. The predicted results for the 
large scale field are compared with the experimental results of Comte-Bellot 
and Corrsln (1971) filtered with (6.1). We find that the averaging scale 
A must be twice the computational mesh scale A for a satisfactory 

a 

calculation of the spectral evolution. We find that calculation In a mesh 

3 

containing as few as 16 points gives remarkably good spectral predictions; 

3 

better results are obtained with 32 points, and It Is reassuring that 
the same constants B, or 6 , fit both sizes. The skewness, which is 

3 

dominated by smaller scales, is predicted much more accurately In the 32 

calculation. Good results are obtained with both (6.3a) and (6.3b). Fig. 4 

3 

shows the results for the 16 calculation. On the basis of this work, we 

now use B^ = 0.06 or B 2 = 0.09 . It Is surprising that 82 ^ 6 ^ , 

because, as Tennekes and Lumley (1972 - Eqn. 3.3.44) show, for 

large Ry . This paradox remains to be understood. 

Next we simulated the TR flow, first with an Initial distribution that 

matched the anisotropy of the TR flow and later with an Isotropic Initial 

distribution. One has problems in setting anisotropic initial conditions 

that are free of shearing stresses, and so the isotropic starting Is probably 

a better approach. It is remarkable that the salient features of the TR 

3 

experiments are captured quite well in a computation using only 16 points! 
The results are shown In Fig. 5; the calculation was executed on a CDC 7600, 
using 120 time steps, In approximately 5 minutes. 

Ferzlqer and Shaanan (1975) are experimenting with a staqqered qrid 
approach that is second-order accurate and does not require explicit 
inclusion of the Leonard stresses. They have validated the constants B 1 
and B 2 with this method, and have also explored the CHC flows. There 
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are some difficulties In providing suitable initial conditions, so comoarison 

with experiments Is not easy. Nevertheless, the salient features of 

3 

the CHC flow can be produced with 16 calculations! 

We have started work on the two-stress mixing layer, in which we expect 
to find a sharp boundary between the turbulent field and irrotational external 
flow. The vorticlty model (6.3b) will offer the advantage of yielding zero 
subgrid scale stresses in a vortici ty-free region, and this is the reason 
that it Is of interest. We plan to extend the computation to "Infinity" 
using inviscid flow theory, and the matching of the computation with the 
inviscid analysis will require use of a difference scheme that does not 
produce vorticlty Improperly. It will be some time before we will feel 
prepared to handle a wall flow accurately. In the meantime, the simulations 
of Deardorff (1970) and Schumann (1973) provide some initial experience with 
channel flows. 

One objective of this work is to test the turbulence models, particularly 
the stress equation model. We can compute the pressure strain terms directly 
(both and T^), and are doing thi« presently. We had hoped that 

the calculations would serve as a basis for evaluating constants in the 
stress-eguation models; instead they seen to be highlighting the weaknesses 
of these models, as discussed in ^5 . However, the fact that a very coarse 
grid produces such remarkably good results leads us to believe that larae 
eddy simulations might, after considerable development, eventually be useful 
for actual engineering analysis. Interested readers should also follow the 
work of Orzag ana Israelii (1972), who are carrying out similar calculations 
usi,.g Fourier ratner than grid methods. 
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